library(dataRetrieval) o2data_raw <- readNWISuv( siteNumbers='02336300', parameterCd='00300', startDate='2013-01-01', endDate='2017-06-18', tz='UTC') o2data <- renameNWISColumns(o2data_raw)
June 19, 2017
Get: dataRetrieval for NWIS/WQP data
Get: dataRetrieval for NWIS/WQP data
## agency_cd site_no dateTime DO_Inst DO_Inst_cd tz_cd ## 1 USGS 02336300 2013-01-03 05:00:00 9.3 A UTC ## 2 USGS 02336300 2013-01-03 05:15:00 9.4 A UTC ## 3 USGS 02336300 2013-01-03 05:30:00 9.4 A UTC ## 4 USGS 02336300 2013-01-03 05:45:00 9.4 A UTC ## 5 USGS 02336300 2013-01-03 06:00:00 9.4 A UTC ## 6 USGS 02336300 2013-01-03 06:15:00 9.5 A UTC ## 7 USGS 02336300 2013-01-03 06:30:00 9.5 A UTC ## 8 USGS 02336300 2013-01-03 06:45:00 9.5 A UTC ## 9 USGS 02336300 2013-01-03 07:00:00 9.5 A UTC ## 10 USGS 02336300 2013-01-03 07:15:00 9.4 A UTC
Get: geoknife for spatial data
library(geoknife)
stencil <- simplegeom(data.frame(c(-84.40764,33.82031)))
daymet <- 'https://thredds.daac.ornl.gov/thredds/dodsC/daymet-v3-agg/na.ncml'
fabric <- webdata(url=daymet, variable='srad', times=c('2013-01-01','2017-01-01'))
job_out <- geoknife(stencil, fabric, REQUIRE_FULL_COVERAGE='false', wait=TRUE)
gdpdata <- result(job_out, with.units = TRUE)
names(gdpdata)[2] <- 'srad'
Get: geoknife for spatial data
## DateTime srad variable statistic units ## 1 2013-01-01 224.0 srad MEAN W/m2 ## 2 2013-01-02 124.8 srad MEAN W/m2 ## 3 2013-01-03 195.2 srad MEAN W/m2 ## 4 2013-01-04 278.4 srad MEAN W/m2 ## 5 2013-01-05 297.6 srad MEAN W/m2 ## 6 2013-01-06 236.8 srad MEAN W/m2 ## 7 2013-01-07 297.6 srad MEAN W/m2 ## 8 2013-01-08 284.8 srad MEAN W/m2 ## 9 2013-01-09 214.4 srad MEAN W/m2 ## 10 2013-01-10 198.4 srad MEAN W/m2
Explore: dygraphs for interaction
library(dygraphs) library(xts) o2data_xts <- xts(o2data$DO_Inst, o2data$dateTime) d <- dygraph(o2data_xts)
Explore: dygraphs for interaction
Explore: ggplot for complexity
library(ggplot2) g <- ggplot(o2data, aes(x=dateTime, y=DO_Inst, color=DO_Inst_cd)) + geom_line() + theme_bw()
Explore: ggplot for complexity
Clean: remove flagged values
library(dplyr) o2data_accepted <- filter(o2data, DO_Inst_cd=='A')
Clean: remove flagged values
## agency_cd site_no dateTime DO_Inst DO_Inst_cd tz_cd ## 1 USGS 02336300 2013-01-03 05:00:00 9.3 A UTC ## 2 USGS 02336300 2013-01-03 05:15:00 9.4 A UTC ## 3 USGS 02336300 2013-01-03 05:30:00 9.4 A UTC ## 4 USGS 02336300 2013-01-03 05:45:00 9.4 A UTC ## 5 USGS 02336300 2013-01-03 06:00:00 9.4 A UTC ## 6 USGS 02336300 2013-01-03 06:15:00 9.5 A UTC ## 7 USGS 02336300 2013-01-03 06:30:00 9.5 A UTC ## 8 USGS 02336300 2013-01-03 06:45:00 9.5 A UTC ## 9 USGS 02336300 2013-01-03 07:00:00 9.5 A UTC ## 10 USGS 02336300 2013-01-03 07:15:00 9.4 A UTC
Analyze: low-O2 events
A common O2 criterion: daily mean >= 5 mg/L, daily minimum >= 4 mg/L
library(dplyr)
o2data_daily <- o2data_accepted %>%
mutate(date=as.Date(dateTime)) %>%
group_by(date) %>%
summarize(DO_Mean=mean(DO_Inst),
DO_Min=min(DO_Inst),
DO_Max=max(DO_Inst)) %>%
mutate(Violation = DO_Mean < 5.0 | DO_Min < 4.0)
Analyze: low-O2 events
g <- ggplot(o2data_daily, aes(x=date)) + geom_ribbon(aes(ymin=DO_Min, ymax=DO_Max), fill='blue', alpha=0.2) + geom_line(aes(y=DO_Mean), color='blue') + geom_point(data=filter(o2data_daily, Violation), aes(y=DO_Mean), color='red') + theme_bw()
Analyze: low-O2 events
Analyze: low-O2 events
o2_violations <- filter(o2data_daily, Violation)
## # A tibble: 17 x 5 ## date DO_Mean DO_Min DO_Max Violation ## <date> <dbl> <dbl> <dbl> <lgl> ## 1 2014-06-21 5.411458 3.9 6.5 TRUE ## 2 2014-06-22 4.927174 3.7 5.9 TRUE ## 3 2014-09-19 5.992708 3.6 8.0 TRUE ## 4 2015-05-26 6.348958 3.1 9.1 TRUE ## 5 2015-05-27 6.955914 3.6 8.4 TRUE ## 6 2015-07-02 5.075789 3.6 6.5 TRUE ## 7 2015-07-23 4.981250 4.2 6.1 TRUE ## 8 2015-07-30 4.954167 4.3 5.8 TRUE ## 9 2015-07-31 4.833333 4.0 5.8 TRUE ## 10 2015-08-02 4.813542 4.3 5.6 TRUE ## 11 2015-08-03 4.733333 4.0 5.6 TRUE ## 12 2015-08-05 5.023656 3.9 5.9 TRUE ## 13 2015-08-06 4.427368 3.5 5.3 TRUE ## 14 2015-08-07 4.925000 4.1 5.4 TRUE ## 15 2015-08-09 4.820000 4.3 5.5 TRUE ## 16 2015-08-10 5.038947 3.8 6.5 TRUE ## 17 2016-10-22 5.426042 3.9 8.1 TRUE
Analyze: streamMetabolizer
nwisdata <- renameNWISColumns(readNWISuv(
siteNumbers='02336300',
parameterCd=c('00300','00060','00010'), # oxygen (mg/L), discharge (cfs), water temperature (C)
startDate='2015-01-01', endDate='2016-01-01', tz='UTC'))
coords <- readNWISsite('02336300')
library(streamMetabolizer)
metabdata <- nwisdata %>%
mutate(
solar.time = convert_UTC_to_solartime(dateTime, longitude=coords$dec_long_va, time.type='mean solar'),
discharge = Flow_Inst * 0.0283168, # cfs to cms
depth = calc_depth(discharge), # uses Raymond et al. 2012 to estimate depth from discharge
DO.sat = calc_DO_sat(Wtemp_Inst, calc_air_pressure(elevation=coords$alt_va * 0.3048)), # serious modeling would use better inputs here
light = calc_light(solar.time, latitude=coords$dec_lat_va, longitude=coords$dec_long_va)) %>%
select(solar.time, DO.obs=DO_Inst, DO.sat, depth, temp.water=Wtemp_Inst, light)
# mm <- metab(specs(mm_name('bayes')), data=metabdata)
Analyze: metabolism at scale (Powell Center synthesis)
library(mda.streams)
mm_bayes <-
get_ts(
c('sitedate_calcLon','gpp_estBest','er_estBest','K600_estBest'),
site_name='nwis_02336300') %>%
unitted::v() %>%
filter(sitedate >= '2013-01-01', sitedate <= '2017-01-01')
Analyze: metabolism at scale (Powell Center synthesis)
library(cowplot)
g <- plot_grid(
ggplot(mm_bayes, aes(x=sitedate, y=gpp)) +
geom_hline(yintercept=0, color='darkgrey') +
geom_point(color='#007929', alpha=0.2) +
theme_bw() + xlab('') + ylab(expression(GPP~(gO[2]~m^{-2}~d^{-1}))),
ggplot(mm_bayes, aes(x=sitedate, y=er)) +
geom_hline(yintercept=0, color='darkgrey') +
geom_point(color='#A64B00', alpha=0.2) +
theme_bw() + xlab('') + ylab(expression(ER~(gO[2]~m^{-2}~d^{-1}))),
nrow=2, align='v')
## Warning: Removed 245 rows containing missing values (geom_point). ## Warning: Removed 245 rows containing missing values (geom_point).
Analyze: metabolism at scale (Powell Center synthesis)
Analyze: rolling means
mm_bayes_rollmean <- mm_bayes %>%
filter(., complete.cases(.)) %>%
mutate(
gpp_roll = zoo::rollmean(gpp, k=15, na.pad=TRUE),
er_roll = zoo::rollmean(er, k=15, na.pad=TRUE),
K600_roll = zoo::rollmean(K600, k=15, na.pad=TRUE))
Analyze: rolling means
g <- plot_grid(
ggplot(mm_bayes_rollmean, aes(x=sitedate, y=gpp)) +
geom_hline(yintercept=0, color='darkgrey') +
geom_point(color='#007929', alpha=0.2) +
geom_line(aes(y=gpp_roll), color='#007929', na.rm=TRUE, size=1) +
theme_bw() + xlab('') + ylab(expression(GPP~(gO[2]~m^{-2}~d^{-1}))),
ggplot(mm_bayes_rollmean, aes(x=sitedate, y=er)) +
geom_hline(yintercept=0, color='darkgrey') +
geom_point(color='#A64B00', alpha=0.2) +
geom_line(aes(y=er_roll), color='#A64B00', na.rm=TRUE, size=1) +
theme_bw() + xlab('') + ylab(expression(ER~(gO[2]~m^{-2}~d^{-1}))),
nrow=2, align='v')
Analyze: rolling means
Analyze: metabolism vs light
ggplot(gdpdata, aes(x=DateTime, y=srad)) + geom_point()
Analyze: metabolism vs light
lightmetab <- inner_join(mutate(gdpdata, date=as.Date(DateTime)), mm_bayes_rollmean, by=c(date='sitedate'))
Analyze: metabolism vs light
## DateTime.x srad variable statistic units date ## 1 2013-01-04 278.4 srad MEAN W/m2 2013-01-04 ## 2 2013-01-05 297.6 srad MEAN W/m2 2013-01-05 ## 3 2013-01-07 297.6 srad MEAN W/m2 2013-01-07 ## 4 2013-01-08 284.8 srad MEAN W/m2 2013-01-08 ## 5 2013-01-09 214.4 srad MEAN W/m2 2013-01-09 ## 6 2013-01-11 121.6 srad MEAN W/m2 2013-01-11 ## 7 2013-01-12 153.6 srad MEAN W/m2 2013-01-12 ## 8 2013-01-13 214.4 srad MEAN W/m2 2013-01-13 ## 9 2013-01-18 240.0 srad MEAN W/m2 2013-01-18 ## 10 2013-01-19 336.0 srad MEAN W/m2 2013-01-19 ## DateTime.y gpp er K600 gpp_roll er_roll ## 1 2013-01-04 17:36:42 -0.06001882 -5.350751 8.031714 NA NA ## 2 2013-01-05 17:36:42 0.10317484 -6.242426 8.494305 NA NA ## 3 2013-01-07 17:36:42 0.19793295 -4.061600 7.016444 NA NA ## 4 2013-01-08 17:36:42 0.09559557 -4.012591 6.702252 NA NA ## 5 2013-01-09 17:36:42 0.14706421 -4.081570 7.311908 NA NA ## 6 2013-01-11 17:36:42 0.35984947 -4.669400 6.816117 NA NA ## 7 2013-01-12 17:36:42 0.38414785 -5.166045 8.194363 NA NA ## 8 2013-01-13 17:36:42 0.47046427 -3.718497 6.200358 0.2430295 -3.962039 ## 9 2013-01-18 17:36:42 0.12745140 -3.583609 5.265111 0.2773623 -3.763658 ## 10 2013-01-19 17:36:42 0.29362743 -4.549904 8.269185 0.2752185 -3.612694 ## K600_roll ## 1 NA ## 2 NA ## 3 NA ## 4 NA ## 5 NA ## 6 NA ## 7 NA ## 8 7.339964 ## 9 7.297291 ## 10 7.161015
Analyze: metabolism vs light
g <- ggplot(lightmetab, aes(x=srad, y=gpp)) +
geom_point(na.rm=TRUE) +
xlab(expression('Shortwave radiation'~(W~m^{-2}))) +
ylab(expression('GPP'~(gO[2]~m^{-2}~d^{-1})))
